#ifndef _BIO_CPP_PARSEFASTA_H
#define _BIO_CPP_PARSEFASTA_H
#include <iostream>
#include <string>
#include <vector>
#include <assert.h>
#include <boost/regex.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/foreach.hpp>
#include "bcpp_parse_error.h"

/*
	Example:

	struct TSequenceData
	{
		list<string>	sequence;
		list<string>	acc_line;
		unsigned		count_seq;
		unsigned		count_invalid;
		bool operator()(const string& accession, const string& sequence, unsigned linecnt);
		bool operator()(const string& accession, const std::vector<std::string>& sequence,
																sequence, unsigned linecnt);
		bool operator()(const string& accession, const string& sequence,
						unsigned acc_filepos, unsigned seq_filepos);
						
	};

	bool process_sequence::operator()(	const string& accession,
										const string& sequence, unsigned line_cnt)
	{
		// reject sequences with a null accession code
		if (acc_line == "SP|000")
		{
			++count_invalid;
			cerr << "Invalid sequence on line "  << line_cnt << std::endl;
			return false;
		}

		++count_seq;
		if (!(count_seq % 100))
			cerr << ".";
		sequence.insert(sequence);
		acc_line.insert(accession);
		return true;
	}

	bool process_sequence::operator()(	const string& accession,
										const string& sequence,
										unsigned acc_filepos,
										unsigned seq_filepos)
	{
		// reject sequences with a null accession code
		if (acc_line == "SP|000")
		{
			++count_invalid;
			cerr	<< "Invalid sequence accession and sequence starting at line: ["
					<< line_cnt << "]"	<< std::endl;
			return false;
		}

		++count_seq;
		if (!(count_seq % 100))
			cerr << ".";
		sequence.insert(sequence);
		acc_line.insert(accession);
		return true;
	}

	bool process_sequence::operator()(	const std::string&				accession,
										const std::vector<std::string>&	sequence,
										unsigned						line_cnt)
	{
		// reject sequences with a null accession code
		if (acc_line == "SP|000")
		{
			++count_invalid;
			cerr	<< "Invalid sequence accession and sequence starting at line: ["
					<< line_cnt << "]"	<< std::endl;
			return false;
		}

		++count_seq;
		if (!(count_seq % 100))
			cerr << ".";
		for (unsigned i = 0; i < sequences.size(); ++i)
			sequence += sequences[i];
		acc_line.insert(accession);
		return true;
	}


	void FooBar_blah_de_blah()
	{
		TSequenceData seq_data;
		Process_FA_Sequences(cin, seq_data);
		Process_FA_Sequences_tellg(cin, seq_data);
		Process_FA_Sequences_vec(cin, seq_data);
	}
*/
namespace BioCPP
{


    struct sequence_formats_alternative_to_fasta
    {
        std::vector<boost::regex> regexes;
        std::vector<boost::shared_ptr<bcpp_parse_format_error> > errors;
        sequence_formats_alternative_to_fasta()
        {
            regexes.push_back(boost::regex("^CLUSTAL", boost::regex::icase));
            errors.push_back(boost::shared_ptr<bcpp_parse_format_error>(new clustal_format_error()));
            regexes.push_back(boost::regex("^\\d+[ \\t]+\\d+", boost::regex::icase));
            errors.push_back(boost::shared_ptr<bcpp_parse_format_error>(new phylip_format_error()));
        }

        void check(const std::string& preamble_str,
                   const std::string& curr_line, unsigned line_cnt)
        {
            for (unsigned i = 0; i < regexes.size(); ++i)
                if (boost::regex_search(curr_line, regexes[i], boost::match_default))
                {
                    std::vector<std::string> preamble;
                    preamble.push_back(preamble_str);
                    errors[i]->set_data(preamble, curr_line, line_cnt);
                    errors[i]->throw_error();
                }
        }
        void check(std::vector<std::string>& preamble,
                   const std::string& curr_line, unsigned line_cnt)
        {
            for (unsigned i = 0; i < regexes.size(); ++i)
                if (boost::regex_search(curr_line, regexes[i], boost::match_default))
                {
                    errors[i]->set_data(preamble, curr_line, line_cnt);
                    errors[i]->throw_error();
                }
        }
    };


//******************************************************************************

//	Process_FA_Sequences

//		reads each sequence and calculates running totals for statistics
//
//******************************************************************************
template <class Ty_>
void Process_FA_Sequences(std::istream& istrm, Ty_& process_sequence)
{
	// These two variables will be non empty AFTER the first sequence is read
	std::string current_sequence;
	std::string current_accession;
	unsigned line_cnt = 0;
	unsigned old_line_cnt = 0;

    bool is_fasta_format = false;

    sequence_formats_alternative_to_fasta alt_formats;

	// get each line in the FASTA file
	std::string line;
	while (std::getline(istrm, line))
	{
		++line_cnt;

		if (!line.length())
			continue;

		// line with sequence data
		if (line [0] != '>')
		{
            if (!is_fasta_format)
                alt_formats.check(current_sequence, line, line_cnt);
			current_sequence += line;
			continue;
		}

        is_fasta_format = true;

		// line with accession number

		// if a previous sequence has been read in successfully, process that			
		if (current_accession.length() && current_sequence.length())
		{
			process_sequence(current_accession, current_sequence, old_line_cnt);
			old_line_cnt = line_cnt;
		}

		// now remember the current accession number and start with an empty sequence			
		// the first character is discarded
		current_accession.assign(line, 1,  std::string::npos);
		current_sequence.erase();

	}


	// if a previous sequence has been read in successfully, process that			
	if (current_accession.length() && current_sequence.length())
		process_sequence(current_accession, current_sequence, old_line_cnt);

	return;
}

template <class Ty_>
void Process_FA_Sequences_vec(std::istream& istrm, Ty_& process_sequence)
{
	// These two variables will be non empty AFTER the first sequence is read
	std::vector<std::string> current_sequence;
	std::string current_accession;
	unsigned line_cnt = 0;
	unsigned old_line_cnt = 0;

    sequence_formats_alternative_to_fasta alt_formats;

    // get each line in the FASTA file
	std::string line;
    bool is_fasta_format = false;
	while (std::getline(istrm, line))
	{
		++line_cnt;

		// line with sequence data
		if (line.empty() || line [0] != '>')
		{
            if (!is_fasta_format)
                alt_formats.check(current_sequence, line, line_cnt);
			current_sequence.push_back(line);
			continue;
		}

        is_fasta_format = true;

		// line with accession number

		// if a previous sequence has been read in successfully, process that			
		if (current_accession.length() && current_sequence.size())
		{
			process_sequence(current_accession, current_sequence, old_line_cnt);
			old_line_cnt = line_cnt;
		}

		// now remember the current accession number and start with an empty sequence			
		// the first character is discarded
		current_accession.assign(line, 1,  std::string::npos);
		current_sequence.clear();

	}


	// if a previous sequence has been read in successfully, process that			
	if (current_accession.length() && current_sequence.size())
		process_sequence(current_accession, current_sequence, old_line_cnt);

	return;
}

//******************************************************************************

//	Process_FA_Sequences_tellg_fast

//		reads each sequence noting positions of accessions and sequences
//
//******************************************************************************
/*
template <class Ty_>
void Process_FA_Sequences_tellg(std::istream& istrm, Ty_& process_sequence)
{
	char buf [32768];
	std::string accession;	accession.reserve(1024);
	std::string sequence;	sequence.reserve(131072);
	bool acc_read = false;
	bool seq_read = false;
	unsigned next_acc_beg = 0;
	unsigned acc_beg = 1;
	unsigned seq_beg = 1;

	char beg_char;
	istrm.read(&beg_char, 1);
	assert(beg_char == '>');
	unsigned buf_pos = 1;


	size_t read_cnt = 0;
	istrm.read(buf, sizeof(buf));
	while ((read_cnt = istrm.gcount()) != 0)
	{
		char* seq_offset = buf;
		char* acc_offset = buf;
		char* buf_end	 = buf + read_cnt;

		// while still can fit both accessions and sequences in buffer
		while (seq_offset && acc_offset)
		{
			//
			// Try to read accession within current buffer
			//
			if (!acc_read && acc_offset)
			{
				// find where sequence begins (after first new line)
				seq_offset = (char*)(memchr(acc_offset, '\n', buf_end - acc_offset));
				
				// not found within this buffer. Append and search in next buffer
				if (!seq_offset)
					accession.append(acc_offset, buf_end - acc_offset);
				
				// found within this buffer. Append. Accession read. Now for sequence..
				else
				{
					accession.append(acc_offset, seq_offset - acc_offset);
					acc_read = true;

					// skip past new line
					++seq_offset;
					seq_beg = (seq_offset - buf) + buf_pos;
					// wrap around
					if (seq_offset == buf_end)
						seq_offset = 0;
				}
			}

			//
			// Try to read sequence within current buffer
			//
			if (acc_read && !seq_read && seq_offset)
			{
				// find where next accession begins (after >)
				acc_offset = (char*)(memchr(seq_offset, '>', buf_end - seq_offset));
				
				// not found within this buffer. Append and search in next buffer
				if (!acc_offset)
					sequence.append(seq_offset, buf_end - seq_offset);
				
				// found within this buffer. Append. Accession read. Now for sequence..
				else
				{
					sequence.append(seq_offset, acc_offset - seq_offset);
					seq_read = true;

					// skip past >
					++acc_offset;
					next_acc_beg = (acc_offset - buf) + buf_pos;
					// wrap around
					if (acc_offset == buf_end)
						acc_offset = 0;
				}
			}

			if (seq_read && acc_read)
			{
				
				assert(sequence.find('>') == std::string::npos);
				assert(accession.find('\n') == std::string::npos);
				if (!process_sequence(accession, sequence, acc_beg, seq_beg))
					return;
				acc_beg = next_acc_beg;
				seq_read = acc_read = false;
				sequence.clear();
				accession.clear();
			}
		}

		// read in next bit of buffer
		buf_pos += read_cnt;
		istrm.read(buf, sizeof(buf));
	}
	
	// if a previous sequence has been read in successfully, process that
	if (acc_read && accession.length())
		process_sequence(accession, sequence, acc_beg, seq_beg);
	

	return;
}
*/





//******************************************************************************

//	Process_FA_Sequences_tellg

//		reads each sequence noting positions of accessions and sequences
//
//******************************************************************************
template <class Ty_>
void Process_FA_Sequences_tellg(std::istream& istrm, Ty_& process_sequence)
{
	char buf [32768];
	std::string accession;	accession.reserve(1024);
	std::string sequence;	sequence.reserve(131072);
	bool acc_read = false;
	bool seq_read = false;
	unsigned next_acc_beg = 0;
	unsigned acc_beg = 1;
	unsigned seq_beg = 1;

    sequence_formats_alternative_to_fasta alt_formats;

    std::string preamble;
    std::string curr_line;
    unsigned line_cnt = 0;
    // skip past everything to first >
	unsigned buf_pos = 0;
	while (istrm && istrm.peek() != '>')
	{
        char c = istrm.get();
        // check each line if if belongs to an alternative format

        if (c == '\n')
        {
            ++line_cnt;
            alt_formats.check(preamble, curr_line, line_cnt);
            preamble += curr_line;
            curr_line = std::string();
        }
        else
            curr_line += c;
		++buf_pos;
	}


	istrm.get();
	++buf_pos;
	// tellg might not work
	//istrm.ignore(std::numeric_limits<int>::max(), '>');

	size_t read_cnt = 0;
	istrm.read(buf, sizeof(buf));
	while ((read_cnt = istrm.gcount()) != 0)
	{
		char* seq_offset = buf;
		char* acc_offset = buf;
		char* buf_end	 = buf + read_cnt;

		// while still can fit both accessions and sequences in buffer
		while (seq_offset && acc_offset)
		{
			//
			// Try to read accession within current buffer
			//
			if (!acc_read && acc_offset)
			{
				// find where sequence begins (after first new line)
				seq_offset = (char*)(memchr(acc_offset, '\n', buf_end - acc_offset));
				
				// not found within this buffer. Append and search in next buffer
				if (!seq_offset)
					accession.append(acc_offset, buf_end - acc_offset);
				
				// found within this buffer. Append. Accession read. Now for sequence..
				else
				{
					accession.append(acc_offset, seq_offset - acc_offset);
					acc_read = true;

					// skip past new line
					++seq_offset;
					seq_beg = (seq_offset - buf) + buf_pos;
					// wrap around
					if (seq_offset == buf_end)
						seq_offset = 0;
				}
			}

			//
			// Try to read sequence within current buffer
			//
			if (acc_read && !seq_read && seq_offset)
			{
				// find where next accession begins (after >)
				acc_offset = (char*)(memchr(seq_offset, '>', buf_end - seq_offset));
				
				// not found within this buffer. Append and search in next buffer
				if (!acc_offset)
					sequence.append(seq_offset, buf_end - seq_offset);
				
				// found within this buffer. Append. Accession read. Now for sequence..
				else
				{
					sequence.append(seq_offset, acc_offset - seq_offset);
					seq_read = true;

					// skip past >
					++acc_offset;
					next_acc_beg = (acc_offset - buf) + buf_pos;
					// wrap around
					if (acc_offset == buf_end)
						acc_offset = 0;
				}
			}

			if (seq_read && acc_read)
			{
				
				assert(sequence.find('>') == std::string::npos);
				assert(accession.find('\n') == std::string::npos);
				if (!process_sequence(accession, sequence,acc_beg,seq_beg))
					return;
				acc_beg = next_acc_beg;
				seq_read = acc_read = false;
				sequence.clear();
				accession.clear();
			}
		}

		// read in next bit of buffer
		buf_pos += read_cnt;
		istrm.read(buf, sizeof(buf));
	}
	
	// if a previous sequence has been read in successfully, process that
	if (acc_read && accession.length())
		process_sequence(accession, sequence,acc_beg,seq_beg);
	

	return;
}

inline void get_seq_at_filepos(std::istream& istrm, unsigned filepos,
							std::string& accession, std::string& sequence)
{
	assert(istrm);
	istrm.seekg(filepos);
	
	// get accession
#ifndef NDEBUG
	istrm.seekg(filepos-1);
	char FASTA_delimiter;
	istrm.get(FASTA_delimiter);
	assert(FASTA_delimiter == '>');
#endif
	getline(istrm, accession);
	
	// get sequence up to next '>'
	getline(istrm, sequence, '>');
	assert(istrm && sequence.length());
}



}

#endif
